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ABSTRACT 

In Smoothed Particles Hydrodynamics (SPH) codes with a large number of parti- 
cles, star formation as well as gas and metal restitution from dying stars can be 
treated statistically. This approach allows to include detailed chemical evolution and 
gas re-ejection with minor computational effort. Here we report on a new statistical 
algorithm for star formation and chemical evolution, especially conceived for SPH 
simulations with large numbers of particles, and for parallel SPH codes. 

For the sake of illustration, we present also two astrophysical simulations obtained 
with this algorithm, implemented into the Tree-SPH code by Lia & Carraro (2000). 
In the first one, wc follow the formation of an individual disc-like galaxy, predict the 
final structure and metallicity evolution, and test resolution effects. 
In the second one we simulate the formation and evolution of a cluster of galaxies, 
to demonstrate the capabilities of the algorithm in investigating the chemo-dynamical 
evolution of galaxies and of the intergalactic medium in a cosmological context. 

Key words: hydrodynamics - methods:numerical - stars: formation - galaxies: evo- 
lution - galaxies: chemical evolution 



1 INTRODUCTION 

The evolution of chemical abundances in the Universe is 
nowadays a subject of utmost importance. Relevant issues 
are when and where the first metals were synthetized, and 
how they spread to enrich the individual host galaxies, 
as well as the intergalactic (IGM) and intracluster (ICM) 
medium; see Ferrara et al. (2000) and Chiosi (2000) for re- 
cent papers on the subject. 

To assess these issues, related to metal production and 
mixing into the interstellar medium (ISM) , it is necessary to 
couple chemical evolution with hydrodynamical evolution. 
This has been done several times in the past, with differ- 
ent techniques. Although in this paper we specifically deal 
with chemical evolution in Smoothed Particle Hydrodynam- 
ics (SPH), we like to recall also some works based on eulerian 
codes for the hydrodynamics. The production and distribu- 
tion of metals over cosmic scales has been recently addressed 
by Cen & Ostriker (1999) and by Yepes et al. (1998), who 
implemented elementary prescriptions for metal enrichment 
in their cosmological simulations. Chemo-dynamical mod- 
els for individual galaxies, considering a multi-phase ISM 
have been developed by Burkert & Hensler (1988), Burk- 
ert et al. (1992), Theis et al. (1992), Samland et al. (1997). 



More recently, Recchi et al. (2001) implemented the detailed 
production and mixing of many different chemical elements 
in their hydro-code for dwarf starburst galaxies. 

The first attempt to couple SPH with chemical evolu- 
tion is by Steinmetz & Miiller (1994), who study the for- 
mation of a Milky Way-like galaxy in a cosmological con- 
text. Their chemical evolution scheme is rather simple: they 
consider the ejection of gas and global metallicity Z from 
Type II supernovae (SN II), occurring over 30 Myr after the 
star formation (SF) episode. These gas and metals produced 
by a star particle are then distributed around over the neigh- 
bouring gas particles. 

More detailed schemes have been worked out by Rai- 
teri et al. (1996) and Berczik (1999). Raiteri et al. model the 
evolution of oxygen and iron, making use of fitting formulas 
to follow gas and metal ejection from a star particle over 
time. Moreover, they include the delayed iron production 
by Type la SN (SN la), by implementing the theoretical 
rate by Greggio & Renzini (1983, hereinafter GR83). The 
gas and metals produced by a star particle are distributed 
around over the neighbouring gas particles, as in Steinmetz 
& Miiller (1994). Berczik follows the method developed by 
Raiteri et al., adding planetary nebulae, hydrogen and he- 
lium to the picture. From the chemical point of view, the 
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main improvement of these models is that they avoid the 
Instantaneous Recycling Approximation (IRA, e.g. Tinsley 
1980) by taking into account the different stellar sources and 
production timescales of the various elements. 

Carraro et al. (1998) and Buonomo et al. (2000) con- 
sider the evolution of the overall metallicity in the IRA, with 
a scheme similar to Steinmetz & Miiller (1994). However, 
in the feedback computation they also include the effect of 
SN la according to the rate by GR83, and they introduce 
a metal diffusion mechanism instead of the standard SPH 
smoothing of metallicity among gas particles. 

Recently, Mosconi et al. (2001) presented a new imple- 
mentation of chemical evolution in SPH, following the evolu- 
tion of very many elements. They consider metal production 
from SN II, essentially instantaneous, and the contribution 
of SN la. The latter, rather than following the theoretical 
rate and time distribution by GR83, is treated simply as 
a prompt metal release occurring after a given time-delay 
tsNi after the SF episode. The typical time-delay tsNi, as 
well as the number of SN la with respect to SN II, are in- 
troduced as free parameters in their simulations. 

All the above mentioned models display one or both of 
the following numerical drawbacks: 

(i) New collisionless star particles are created when 
SF occurs, while the corresponding metal ejecta are redis- 
tributed around to other gas particles. Hence, the number 
of SF episodes and of "offspring stars" per gas particle must 
be artificially damped, otherwise the overall number of par- 
ticles gets too large. Also, the simulation must start with 
relatively few baryonic particles, as SF will increase their 
number substantially in the course of the simulation; this 
hampers the modelling of the early, very interesting phases 
of galaxy formation. 

(ii) When SF occurs, hybrid particles are created which 
host both gas and different stellar populations. As hybrid 
particles consist of both collisional and collisionless sub- 
components, from the point of view of the hydrodynamics 
they introduce artifacts: before becoming purely collisionless 
particles, the stars tend to follow for a while the dynamical 
evolution of the gas, clearly a spurious effect. 

In this paper, we present a new statistical algorithm 
of SF and chemical evolution in the context of N-body SPH 
simulations, aimed at overcoming the above mentioned prob- 
lems. In our approach, in fact, the number of baryonic par- 
ticles remains constant throughout the simulation, and par- 
ticles are either fully hydrodynamical or fully collisionless. 

Due to its low computational costs, the new formalism 
is particularly suited to many-particle simulations, and the 
completely local character of the computation of "chemical 
quantities" makes it convenient for parallel codes (Lia 2000, 
Lia & Carraro 2000, Dalla Vecchia 2001, Lia et al. 2001). 

The scheme follows the detailed production of many 
chemical elements over different timescales, avoiding the 
IRA. In this respect, we provide updated fitting formula? for 
the evolution of gas and metal release, that can be easily im- 
plemented also in other particle-based codes. The approach 
allows to calculate the abundance evolution of a number of 
chemical elements independently, as well as of the global 
metallicity. Depending on the objects being modelled, one 
might in fact be interested in monitoring different chemical 
elements in different cases. 

We remark that relaxing the IRA is fundamental not 



only for the sake of chemical evolution (to trace abun- 
dance ratios of different elements, for instance [a/Fe] ra- 
tios), but also for the effects that delayed gas restitution 
may have on the dynamics, as recently discussed by Jung- 
wiert et al. (2001). Among others, continuous mass loss from 
stars over a Hubble time contributes to solve Robert's time 
paradox in spiral galaxies (Kennicutt et al. 1994). 

The layout of the paper is as follows. In Section 2 
we briefly introduce our TreeSPH code, while in Section 3 
we describe the star formation algorithm and related feed- 
back prescriptions. Section 4 provides a detailed explana- 
tion of our statistical chemical model, whereas Sections 5 
to 8 present tests and astrophysical applications. Finally, 
Section 9 summarizes our findings. 



2 THE N-BODY TREE-SPH CODE 

The simulations presented here have been performed using 
the Tree-SPH code developed in Carraro et al. (1998), Lia 
(2000), Lia & Carraro (2000) and Dalla Vecchia (2001). In 
this code, the properties of the gaseous component are de- 
scribed by means of the SPH technique (Lucy 1977; Gin- 
gold & Monaghan 1977), whereas the gravitational forces 
are computed by means of the hierarchical tree algorithm of 
Barnes & Hut (1986), using the tolerance parameter 6 = 0.8, 
expanding the tree nodes to quadrupole order, and adopt- 
ing a spline softening parameter. In the SPH method, each 
particle represents a fluid element whose position, velocity, 
energy, density etc. are followed in time and space. The prop- 
erties of the fluid are locally estimated by an interpolation 
which involves the smoothing length h. Convergence tests 
for the SPH part of the code were performed in Carraro 
et al. (1998). Finally, cooling tables as a function of the 
metallicity of the gas have been taken from Sutherland & 
Dopita (1993). For further details, we refer the reader to the 
above mentioned papers. 

In this paper we improve upon the SF and Chemical 
Evolution sections of our code, as presented here below. 



3 STAR FORMATION AND FEED-BACK 

The formation of stars is a poorly understood process, and 
therefore difficult to model properly. The most widely used 
strategy consists of two steps. 

Firstly, an element of fluid must satisfy some conditions 
to be eligible to SF. The basic criteria adopted in our SPH 
code to select the fluid elements prone to SF are (i) the 
gas particle must be in a convergent flow, and (ii) the gas 
particle must be Jeans unstable. The conditions are met if 
the velocity divergence is negative and if the sound crossing 
time-scale is shorter than the dynamical time-scale. The ef- 
fect of varying the SF criteria has been discussed in details 
in Buonomo et al. (2000), which the reader is referred to. 

Secondly, the fluid element turns into stars according to 
a suitable star formation rate (SFR), which is customarily 
a reminiscence of the Schmidt (1959) law: 

dp* _ _°}Pg_ _ c Pg 
dt dt tff 

where c* is the star efficency parameter, that we set equal 
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to 0.1, and i// is the free-fall timescale. In our code, this 
SF law is seen as the probability that a single gas particle is 
completely transformed into a star particle in the next time 
step (Katz & Gunn 1991; Katz et al. 1996). This probability 
is computed as 



P(SF) = 1 - exp{- 



c*At. 



(1) 



where At is the particle time step. The particle is then ef- 
fectively chosen, or not, to transform into a star particle by 
means of a Monte Carlo method. This statistical approach 
is the more reliable the larger is the number of particles that 
model the star forming region. 

The advantage of this approach is twofold. First of all, 
the total number of particles in the simulation does not 
increase, since a gas particle turns completely into a star 
particle once it fulfills all the SF criteria. Besides, baryonic 
particles are assigned and always retain all the same mass 
throughout the simulation. Secondly, at increasing number 
of particles enclosed in a star forming region, the probability 
approach naturally becomes more and more realistic. Again, 
this makes it potentially advantageous for algorithms run- 
ning on parallel computers, where large numbers of particles 
can be managed (Lia 2000; Dalla Vecchia 2001). 

Stars are expected to return energy to the ISM via 
SN explosions, and the bulk of the energy released in the 
ISM is known as the stellar energy Feed-Back. In our code, 
the SN rates are calculated by the statistical algorithm as 
explained later in this paper, and each SN is assumed to 
contribute 10 50 erg of energy to the surrounding medium 
(Thornton et al 1998). 

The key problem here is to know how the released en- 
ergy is given to the ISM, because the limited resolution of 
N-body simulations does not allow to describe the ISM as a 
multi-phase medium. In Buonomo et al. (2000) we have in- 
vestigated different prescriptions for the Energy Feed-Back, 
suggesting that a good solution is to release all the energy 
from SN and other sources to the thermal budget of the 
fluid element (Steinmetz & Miiller 1994). As in this scheme 
the feed-back is localized within the gas particle, and gets 
rapidly dissipated by cooling, it turns out to have no signifi- 
cant effects on the dynamics and/or on the chemical proper- 
ties of the gas throu gh th e mixing of metals (which is treated 
indipendently, see § 4.4). 



4 THE "STATISTICAL" CHEMICAL 
ALGORITHM 

Once stars are present, they return to the ISM part of their 
mass, in form of chemically processed gas, and energy, via 
the stellar feed-back. In this paper we extend the probabilis- 
tic approach adopted for the SF process to the energy and 
chemical feed-back. Each baryonic particle can be either in 
the form of gas or in the form of stars. The probability that 
a particle turns from gas to stars is fixed by the SFR (Eq. 
Likewise, each star particle can be assigned a certain prob- 
ability to turn back into gas; and when this takes place, the 
particle carries along the corresponding metal production, 
SN rate, and energy feed-back. Notice that in the present 
approach all baryonic particles, stars and gas, have the same 



mass and their total number is conserved. Our statistical al- 
gorithm is developed and only valid under this condition. 

In brief, we consider a star particle as a Single Stellar 
Population (SSP) of assigned Initial Mass Function (IMF), 
and calculate the fraction of its mass which turns into re- 
ejected g SiS . £LS a function of its age. This mass fraction is 
identified as the probability that a star particle of age t be- 
comes a gas particle again at time t + At, being At the 
particle time step. For a Salpeter IMF, for instance, this 
probability is of the order of 30% over a Hubble time, since 
this is the fraction of the initial mass of the SSP which is 
globally re-ejected. By means of a Monte Carlo method, at 
each timestep we transform back into gas a certain frac- 
tion of the star particles; typically a 30% of the total stars 
formed will have become gas again by the end of the sim- 
ulation, after a Hubble time. Over a large enough number 
of particles, star formation and the corresponding gas resti- 
tution should be well represented statistically. Besides, this 
approach easily takes into account that gas restitution from 
stars is diluted in time, with no need to assume the IRA. 

Our statistical chemical algorithm is described in detail 
in the following sections. 



4.1 Gas restitution 

Each star particle is treated as a SSP of total mass m and age 
t = T — To, where T is the present age of the system and To is 
the birth-time of the star particle. Ideally, one could follow 
the detailed chemical, photometric and spectral evolution of 
each SSP by adopting a grid of stellar models (e.g. Chiosi 
et al. 1998). However, in a hydrodynamical code this would 
heavily increase the computational load . To circumvent this, 
we prefer to approximate the amount of released gas, the 
SN rates and the metal production by means of analytical 
formulas, as described here below. 

Within a SSP, stellar masses are distributed ac- 
cording to the IMF $(M), usually a power-law (e.g. 
$(M) cx M -1,35 for a Salpeter IMF). A star of given mass M 
is characterized by a lifetime r(M); when it "dies", part of 
it remains enclosed in a remnant of mass M r (a white dwarf, 
a neutron star or a black hole) while the rest is ejected back 
to the ISM in the form of chemically processed gas. Stellar 
models provide r and M r as a function of M. The mass 
fraction of the SSP which is ejected back in the form of gas 
by time T = T + 1 is: 



E(t) = 



where 



e(t) = 



A/„. 



M -M r (M) 



M(t) 



M 



$(M)dM= / e(t')dt 



M — M r (M) 
M 



dM\ 



(2) 



M(t) 



M u — 100 Mq is the upper stellar mass limit in the IMF, 
M(t) indicates the mass of a star with lifetime r, and 
M(t) = M(T - T ) is thus the smallest stellar mass dy- 
ing and restituting gas by time T (see also Tinsely 1980). 
For a Salpeter IMF, for instance, the global returned gas 
fraction E over a Hubble time is of the order of 30%. 

Accordingly, e(t) of Eq. (Q) is the rate of gas ejection in 
time. In our statistical approach, this rate is interpreted as 
the probability that a star particle, in time, transforms again 
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Figure 1. Rate of gas restitution for SSPs with different IMFs, 
as a function of the age t of the SSP. Thick lines: numerical re- 
sults; thin lines: corresponding power-law fits. The function e(t) 
represents as well the probability per year that a star particle 
transforms back into a gas particle (see text). 



into gas. Namely, at each time step At the star particle is 
assigned a probability 

/t+Ai 
e(t') dt' (3) 

to turn back into gas, and a Monte Carlo method decides 
whether it does switch to a gas particle within At. The over- 
all probability that a star particle ever switches to gas again, 
integrated over time, is E ~ 30%; accordingly, of all the star 
particles formed at time To roughly one third will have re- 
turned to be gas after a Hubble time. Over a sufficiently 
large number of particles, this approach should give a fair 
representation not only of the overall gas restitution, but 
also of its rate in time (avoiding the IRA). 

Actually if, say, N particles become star particles at 
some time To, to yield the correct gas restitution in time the 
probability (H) should be applied to all of those N particles, 
at each timestep, throughout the simulation. However, by 
age t a fraction E(t) of those N particles will have already 
returned to be gas, while only a fraction 1 — E(t) will still be 
stars. Obviously, the probability to transform back into gas 
at age t can be calculated only for those N [1— E(t)] particles 
that are still stars at that time, rather than on the base of the 
whole initial population of N particles. Hence, to recover the 
correct statistical gas restitution, the probability (|^) must 
be corrected by a further factor [1 — E(t)\, and becomes: 

J t t+At e(t')dt' 



gt (At) = 



(4) 



l-E(t) 

Our test in § | and Fig. g demonstrate that in this way 
the correct global gas restitution is statistically recovered 
- which would not be the case if expression Q for the 
probability were directly used. 

We should now determine e(t) as a function of the age t 
of the SSP. As mentioned above, in terms of computational 
effort it is convenient to adopt some analytical fit to the 
numerical results of detailed chemical evolution models. We 
proceed as follows. 
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Table 1. Rate of gas restitution e(t) for SSPs with different IMFs; 
e(t) = for t < 3.4 Myr, lifetime of the most massive star, 
M u = 100 Mq. e(t) represents as well the probability per year 
that a star particle transforms back into a gas particle (see text). 



We calculate numerically the function adopting the 
stellar lifetimes r(M) and remnant masses M r (M) from the 
Padova evolutionary tracks (see Portinari et al. 1998, here- 
inafter PCB98, and references therein). Regarding the IMF, 
its shape and variations with the environment are still an 
open issue (e.g. Scalo 1998); hence we consider three differ- 
ent IMFs suggested in literature, with the intent of covering 
a representative range of possibilities. 

(i) The standard and widely adopted Salpeter (1955) IMF 
$(M) = C s Af -135 C s = 0.1716 

(ii) The somewhat steeper IMF by Kroupa (1998), best 
suited to the Solar Neighbourhood: 

C C kl M~ - 5 
#(M)= I C ft M- 12 
I C k M- 17 



M < 0.5 
0.5 < M < 1 



C k i = 0.48 



M > 1 
C k = 0.295 



(iii) The more top-heavy IMF suggested by Arimoto & 
Yoshii (1987) for elliptical galaxies: 



$(M) = C a M~ 



C a = 0.145 



C s , C a , Ck and Cki are normalization coefficients fixed so 
that the IMF is normalized to unit mass when integrated 
between the low and high stellar mass ends (0.1 and 100 Mq 
respectively). With the adopted stellar mass limits, SSPs 
with the Salpeter or Kroupa IMF restitute a 30% of their 
initial mass in gas, while for SSPs with the more top-heavy 
IMF by Arimoto & Yoshii the restitution fraction is larger 
(about 50%). 

The numerical results for the function (^) are well fit- 
ted by power laws (Fig. ^|) ; the fitting functions are listed in 
Table m In principle, stellar lifetimes and remnant masses 
depend also on the metallicity of the star (PCB98), and 
hence e(t) also depends on the metallicity of the parent SSP. 
However, it is beyond the scope of this paper to implement 
chemical evolution in such fine detail; therefore, for gas resti- 
tution (and for the chemical yields below) we adopt fits in- 
dependent of metallicity. Notice, however, that metallicity- 
dependent prescriptions can in principle be inserted in our 
approach, as a star particle always keeps track of the metal- 
licity Zo of its host SSP (as we did for the sole case of the 
nitrogen yields in 



4.3.1 



With the analytical expressions given for e(t) in Table |l|, 
once the IMF is chosen it is straightforward to calculate, 
at each time step, the probability (^) that a star particle 
transforms back into gas. If this happens, in the following 
discussion we will call it a "gas-again particle" . 
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Table 2. SN rates for SSPs with different IMFs; the rates of SN II and of SN la drop to zero out of the respective time ranges of activity, 
indicated in the rightmost column. 
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Figure 2. Rate of SN II for SSPs with different IMFs, as a func- 
tion of the age t of the SSP. Thick lines: numerical results; thin 
lines: corresponding power— law fits. 



4.2 The supernova rate 

The SN rate enters the calculation of chemical evolution and 
of feed-back. We determine the SN rate considering each 
particle that undergoes SF as a SSP. We discuss the rates 
of supernova? type II and type la separately. 



4-2.1 The rate of SN II 

SN II (and SN Ib,c) originate from single stars of mass 
M > M up , the minimum stellar mass allowing for non- 
degenerate carbon ignition in the core; we refer to the 
Padova stellar tracks, where M up = 6 M . As $(M)/M is 
the distribution of stellar masses by number of stars, a SSP 
of unit total mass produces a rate of SN II in time given by: 



m(t) = 



M(t) 



r (M u ) < t < r(M up ) 



otherwise 



Just as we did for the gas restitution fraction, we calculate 
rn(t) numerically for our three IMFs and fit the numer- 
ical results with suitable analytical functions (Fig. |^ and 
Table |). 



A particle of mass m which underwent SF produces in 
a time step At a number of SN II given by: 

/t+At 
rii(t')dt' (5) 

The factor m must be introduced because rn(t) is the SN 
rate for a SSP of total unit mass (from the normalization 
of the IMF), while our baryonic particles in general have a 
total mass m. 



4.2.2 The rate of SN la 

For the rate of SN la we follow the scheme of GR83, assum- 
ing SN la to originate from binary systems of total mass Mb 
between 



M BJnf = 3M( 







and 



Mb.bup = 2 M up = 12M, 







over timescales set by the lifetime of the secondary star M2 
in the system. In this scenario, the total number of SN la 
produced by a SSP of age t is: 



RsNIa(t) = A 



where 



/(/*) = 24 M 2 



_ Mb~ 



(1Mb 



Mo 

M=^€ [0,0.5] 



Hi„f — max 



M 2 (t) M B -M u 
Mb ' Mb 



(see GR83, for further details). A is a parameter fixed so as 
to match the observed ratio of SN II/SN la in Milky Way- 
type galaxies; for M B , in f = 3 M , A ~ 0.07 (e.g. PCB98). 

Since the timescale of explosion is set by the lifetime of 
the secondary star, to calculate the rate of SN la in time we 
first invert the order of integration: 



RsNIa{t) = A 



M 2 ,$(Mb] 



Me 



Ml 



dM E 



dM 2 



'M 2 {t) 

where 

M B ,m,in = max{M B ,i„/, 2M 2 } 
M B , m ax = min{M fl , s „p, M 2 + M up } 

Notice that the formula has been corrected by a factor Mg 1 
in the inner integral, with respect to the analogous formula 
given for a single SF burst by GR83. 



© 2001 RAS, MNRAS 000, 000-000 



6 Lia, Portinari and Carraro 



Salpeter 
Arimoto — Yoshii 
Kroupa 




Figure 3. Rate of SN la for SSPs of age t and with different 
IMFs. 



Accordingly, the rate of SN la in time is: 



dM 2 \ f M B, 



/( 



M 2 s HMb) 



ria(t) = 



dM B 

J M 2 (t) 

t > r(M up ) 
otherwise 



We calculate rj a (t) numerically for our three IMFs and pro- 
vide analytical fits to the numerical results (Fig. ^ and Ta- 
ble |). 

A particle of mass m which underwent SF produces in 
a time step At a number of SN la given by: 



ft+Ai 

N S Nia = m I n a (t')dt' 



4-2.3 The statistical correction 



(6) 



So far, we have derived the SN rates for a SSP, and conse- 
quently for a particle which has experienced SF. When we 
are to calculate the feed-back effect that supernovae have on 
a given gas particle, in principle we need to know the number 
of SN produced in the time step At by all of its neighbouring 
particles which have experienced SF at some previous time. 
These include both star particles and gas-again particles. 
Although for each of these particles it is straightforward to 
calculate the corresponding Nsnii @ and Nsnio. (|^), the 
above mentioned procedure would require to determine all 
the neighbouring particles, both gas and stars. On the other 
hand, for hydrodynamical purposes only the neighbouring 
collisional (gas and gas-again) particles are of interest. So, in 
principle a double calculation of neighbours would be neces- 
sary, which would require a major computational effort and 
could not be trivially implemented in parallel codes. 

We thus restrict the calculation of neighbours to gas 
particles, as usual in SPH codes. As a consequence, we can 
determine the contribution to the SN rate (and to the related 
feed-back, see Section 4) only from neighbouring gas-again 
particles; we should therefore correct for the "missing" con- 
tribution of star particles. We proceed as follows. 



SN rates within the timestep At are calculated, rather 
than for all the particles that have experienced SF, only for 
those particles that return to be gas right within At. These 
particles represent a fraction gt(At) — given by (^) — of the 
whole parent stellar population. Accordingly, the "statistical 
correction" must be 



N S N 



N SN 
9t(At) 



(7) 



Notice that in this approach not all gas-again particles 
(namely, any gas particles that have ever experienced SF in 
their past) contribute to the SN rates, but only those that 
become gas-again "right now" . This ensur es tha t the SN ex- 
plosions (and the metal production, see 



4.3.1 



below) take 

place exactly where the parent stellar component is located. 

The statistical correction (^) is quite large: a whole stel- 
lar population is sampled by the fraction of it which dies 
within the timestep At. This fraction can become very small, 
especially for small timesteps and/or at advanced ages after 
the SF episode (see the behaviour of e(t) in time, Fig. Q). 
At late times, the hydrodynamical timestep of the system is 
very small with respect to the rate of gas restitution, hence 
statistical fluctuations in the number of old star particles 
turning to gas-again tend to become quite large. This might 
induce large fluctuations in the effective release of SN and 
metals with respect the smooth theoretical time evolution 
of a SSP. To hamper these fluctuation, the probability that 
a star particle becomes gas-again is computed, rather than 
over the dynamical timestep, over a "chemical timestep" 
which increases in proportion with the age of the SSP, so 
as to smoothen the related increase of the statistical noise. 
For the sake of clarity, from here on we indicate by At and 
St the chemical and dynamical timestep respectively. The 
chemical timestep At is chosen to be the minimum between 
St and one tenth of the age t of the SSP. We also put an 
upper limit to At of 2 x 10 s yr, so that it does not become 
excessively large from the dynamical point of view, when we 
are dealing with stars that are a few Gyr old. 2 x 10 8 yr 
corresponds, for instance, to the typical revolution period of 
the Sun in the Galaxy. 

In detail, our algorithm works as follows. For a star par- 
ticle, the "Monte Carlo chance" to transform to gas-again 
is activated once every At, rather than at each dynamical 
step St, with a corresponding probability (^) also calculated 
for At. If the star particle is selected to become a gas-again 
particle, it "produces" a number of SN (fjj) also calculated 
on the base of At. Before the next dynamical timestep, this 
particle is added to the list of gas particles for the calcula- 
tion of neighbours, and from then on it behaves just like any 
other gas particle. Finally, it is worth noting that the energy 
feed-back is released over the hydrodynamical timestep, and 
the energy balance between cooling and heating is calculated 
over the hydrodynamical timestep, as well. 

The various test applications presented later in this pa- 
per show that this statistical approach, with the adopted 
chemical timestep, yields sensible results. 



4.3 Chemical enrichment 

The calculation of the chemical enrichment of the gas pro- 
ceeds in two steps: 
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(i) calculation of metal production by SSPs in particles 
which have experienced SF; 

(ii) metal diffusion among gas particles. 

In the tables we provide the necessary information to calcu- 
late the production of various chemical elements, as well as 
that of the overall metallicity, independently of one another. 
In fact, according to the particular object being modelled, 
one might be interested in tracking different chemical ele- 
ments because available abundance data also depend on the 
class of object considered. Hence, in each simulation one can 
select and "switch on" the elements of interest. The hydro- 
gen abundance, necessary to express chemical abundances 
in the usual [element/H] dex scale, can be obtained for each 
particle as: 

X =1-Y-Z 

where Y represent the helium mass fraction and Z the over- 
all metallicity. 

For the sake of clarity, here below we will describe our 
implementation of chemical evolution in terms of a single 
"metal parameter" Z. This parameter is meant to represent 
the chemical abundance of any specific chemical element — 
helium or metals. 



4-3.1 Metal production of a star particle 

Let's consider a gas particle of metallicity Zq which becomes 
a star particle at time To: it effectively hosts a SSP composed 
of a distribution of stellar masses <E>(M), all born at time 
To out of gas with homogeneous metallicity Zq. As time 
progresses, stars of smaller and smaller mass M die, each 
ejecting a mass of metals given by: 

M z =yz + Z (M - M r ) 

The first term yz is the stellar yield, i.e. the newly synthe- 
sized metals (see the definition by Tinsley 1980 as revised 
by Maeder 1992); the second term is the metals present in 
the star from birth and re-ejected. The mass fraction that a 
SSP releases in the form of metals up to age t is: 



E z (t) =1 



M(t) M 



$(M) dM 



=i M ;i $( M)dMiM;;/' f 



r *(m) 



dM ' 
dr , 



t) 
dt' 



£$(M) dM 



Z E(t) 



r(M„) M \ dr } \ M (t>) 

Evidently, summing the Ez's over all the chemical elements 
one expects: 

Y,Ez(t) = E(t) 

The rate of metal ejection becomes: 
ez(t) =pz(t) + Z e(t) 
which is determined once we know 
$(M) / <IM ' 



Pz(t) = 



yz ■ 



M 



/ dM \ 
V "dr J 



(8) 



M(t) 



for the chemical ele ment Z of interest, since e(t) has already 



been calculated in 6 4.1 



To calculate pz(t), we adopt the stellar yields yz by 
PCB98 for the case of massive stars, and by Marigo (2001, 



her case a = 1.68) for low and intermediate mass stars. 
Similarly to what we did for the returned gas fraction and 
for the SN rates, we calculate the functions (^j) numerically, 
for several chemical elements and for the three IMF cases. 
Then we provide suitable analytical fits to the numerical 
results, to be used in the hydrodynamical code. Our fitting 
functions are listed in Table |j. As in the case of the returned 
gas fraction, we neglect in general the dependence of stellar 
yields on the initial metallicity of the SSP. 

The derived fitting functions reflect the different nu- 
cleosynthetic history of the various elements. Both helium 
and metals (meant as the overall metallicity) are expelled 
over the whole mass range; although the bulk comes from 
massive stars, the contribution from stars of intermediate 
and low mass cannot be neglected. Therefore, the "onset" 
of the fitting functions for He and Z corresponds to the life- 
time of the most massive star (3.4 Myr) and the production 
continues forever; calculations have though been stopped at 
15 Gyr, beyond that (which is beyond the age of the Uni- 
verse anyways) the production can be considered negligi- 
ble. The production rates rapidly decrease in time, due to 
the behaviour of stellar lifetimes with mass (cf. the factor 
dM/dr in Eq. ||) . The Arimoto-Yoshii IMF case is the most 
skewed toward massive stars, hence the metal production is 
the largest and it presents the steepest drop with time as it 
is the most dominated by massive stars. The global metal 
production is lower and its decrease with time is slower going 
to the Salpeter and then to the Kroupa IMF case. 

Oxygen production, due to massive stars, is limited to 
the range of lifetimes of SN II progenitors (3.4 to 34 Myr). 
Especially for the most massive stars, oxygen production 
is sensitive to metallicity (PCB98). In the framework of a 
chemical algorithm handy enough to be implemented in hy- 
drodynamical codes, it is not worth entering so much detail 
and we rather give average estimates of oxygen yields. Again, 
the global oxygen production is largest for the Arimoto- 
Yoshii IMF case and decreases going to the Salpeter and to 
the Kroupa case. 

Similar trends, and comments, apply to the production 
of iron, magnesium, silicon, sulfur and calcium, also ejected 
by massive stars. Notice that we are discussing here only 
chemical yields from single stars; the contribution of SN la 
(binaries) will be added later. Theoretical magnesium yields 
are known to be underestimated with respect to observations 
(Timmes et al. 1995, Thomas et al. 1998, PCB98, Chiappini 
et al. 1999) , while iron yields are sometimes suggested to be 
a bit high (though still within the intrinsic uncertainty of 
a factor of 2 in SN models, Timmes et al. 1995). The latter 
point is true especially when considering the latest estimates 
of the solar abundance JFe/H]0=7. 5, rather than the old 
value of 7.67. In Table H, entries for magnesium and iron 
have been optimized on the base of observational indications 
for the Solar Vicinity. 

The nucleosynthetic history of carbon is quite com- 
posite. Carbon is ejected by massive stars through stellar 
winds and SN II explosions, while intermediate and low- 
mass stars contribute to its production in their TP-AGB 
phase. In both cases, theoretical carbon yields are sen- 
sitive to metallicity effects (PCB98, Marigo 2001). From 
the point of view of observations, the role of the vari- 
ous contributors and of metallicity dependence is still de- 
bated (Prantzos et al. 1994, Gustafsson et al. 1999, Gar- 
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chemical element 


Salpeter 


Kroupa 


Arimoto-Yoshii 


time range [yr] 


Zj 


Pz(t) — 


453 1 -1 - 7 


lO.Si" 1 ' 5 


5919 t -1 - 8 


t e [3.410 6 ,1510 9 ] 


4 He 


PHe(t) = 


2.8t -1 - 37 


0.14 1 -121 


80i -152 


f f [3 4 1 6 1 5 1 9 1 


16 Q 


Po(t) = I 


7.3e-10 
7.3e4t" 2 


4.2e-10 
4.2e4t~ 2 


1.8e-9 
1.8e5i~ 2 


t G [3.4 10 6 , 10 7 ] 
t e [10 7 ,3.410 7 ] 


56 Fe 


PFe{t) = < 


5.1e-ll 
5.1e3t~ 2 


3.3e-ll 
3.3e3t~ 2 


l.le-10 
l.le4t" 2 


t e [3.4 10 6 , 10 7 ] 
t e [10 7 ,3.410 7 ] 


24 Mg 


PMg(t) = < 


4.6e-ll 
4.6e3t -2 


2.5e-ll 
2.5e3t -2 


1.2e-10 
1.2e4t" 2 


t e [3.4 10 6 , 10 7 ] 
t e [10 7 ,3.410 7 ] 


28 Si 


PSi(i) = I 


5.8e-ll 
5.8e3i~ 2 


3.6e-ll 
3.6e3i" 2 


1.3e-10 
1.3e4t" 2 


t G [3.4 10 6 , 10 7 ] 
t G [10 7 ,3.410 7 ] 


32g 


Ps(t) = I 


2.9e-ll 
2.9e3i" 2 


1.8e-ll 
1.8e3i~ 2 


6.6e-ll 
6.6e3t" 2 


t G [3.4 10 6 , 10 7 ] 
t G [10 7 ,3.410 7 ] 


40 Ca 


PCa(t) = < 


4.2e-12 
4.2e2i" 2 


2.7e-12 
2.7e2i" 2 


9.1e-12 
9.1e2t" 2 


t G [3.4 10 6 , 10 7 ] 
t G [10 7 ,3.410 7 ] 


12 C 


Pc(t) = } 


1.2e-10 
1.2e4t" 2 
l.le-6t-°' 7 - 3e-13 


5.8e-ll 
5.8e3t~ 2 
2.9e-8t-°' 5 - 8e-13 


3.6e-10 
3.6e4i" 2 
6.4e-5t-°- 9 - 2e-13 


t G [3.410 6 ,10 7 ] 
t G [10 7 ,3.410 7 ] 
t G [210 8 ,510 9 ] 


14 N 


PNs(t) = 


7.7 Z t~ 1A 


0.98 Zo*" 1 ' 3 


520 Zo*" 1 ' 6 


t G [3.410 6 ,1510 9 ] 


PNp(t) = < 


Zo < 0.004 

Z G [0.004, 0.02] 

Z > 0.02 


3.3e-12 

3.3e-18 
Z 

5.8e-14 


3.3e-12 

3.3e-18 
Z 

5.8e-14 


4.6e-12 
4.7e-18 

8.3e-14 


t G [10 8 , 2.5 10 s ] 



Table 3. Rate of release of chemical yields of various elements for SSPs of age t and with different IMFs; the rates drop to zero out 
of the respective time ranges, indicated in the rightmost column. In the case of nitrogen, we distinguish the secondary and the primary 
components and give metallicity-dependent prescriptions; Zo is the initial metallicity of the SSP. Needless to say, the total yield for 
nitrogen is pjv(t) = pjv s (t) + PjVp(i) ( see text for details). 



nett et al. 1999, Henry et al. 2000, Carigi 2000, Liang et al. 
2001). Considering these uncertanties, for the purpose of the 
present chemical algorithm we simply provide carbon yields 
averaged over metallicity both for massive stars (time range 
3.4-34 Myr) and for the delayed contribution of lower mass 
stars (0.2-5 Gyr). Also for carbon the entries in Table ^ 
have been adjusted on the base of observational constraints 
for the Solar Vicinity. 

The stellar nucleosynthesis of nitrogen is also complex, 
and the distinction between its primary and secondary pro- 
duction (Tinsley 1980) is of prime importance to interpret 
observational evidence (e.g. Larsen et al. 2001 and refer- 
ences therein). For nitrogen, therefore, we consider it worth 
treating the secondary and primary components separately 
(pNs and Pn p ), with metallicity dependent prescriptions. 
The secondary component is, by definition, directly propor- 
tional to the metallicity of the parent SSP and is produced 
by stars of all masses (from 3.4 Myr to 15 Gyr). Primary 
production occurs in intermediate mass stars of 3.5-5 Mq 
(time range ~ 100-250 Myr) and is very sensitive to metallic- 
ity, being more efficient at lower metallicities (Marigo 2001, 
her case a=1.68 adopted here). Hence we give a metallicity 
dependent fit also for the primary component. For carbon 
and nitrogen, the relative importance of the production by 
intermediate and low mass stars with respect to that by mas- 



sive stars increases when moving from the Arimoto-Yoshii 
to the Salpeter to the Kroupa IMF, due to the corresponding 
shift of the SSP towards less massive stars. 

The contribution of SN la 

Table ^ provides the rate of metal release in time from sin- 
gle stars in a SSP of total unit mass. We should then add 
the contribution of SN la, which originate in binary sys- 
tems. This is easily done since we know the rate of SN la in 
time (Table ||); for each SN la, we adopt the chemical ejecta 
M' z a from the W7 model by Iwamoto et al. (1999). The to- 
tal metal release by the SSP, including the contribution of 
SN la, is given by: 

p t z° t (t)=Pz(t)+M I z a xr Ia (t) 

as listed in Table ^[ 

Hence, within a timestep At a SSP of age t releases an 
amount of metals given by: 

/t+At fi+At ft+At 

e z (t')dt' = I p^it^dt' + Z / e(t')dt' 

or, equivalently, ejects an amount J t ' +A * e(f') dt' of gas with 
metallicity: 



© 2001 RAS, MNRAS 000, 000-000 



Star formation and chemical evolution in SPH simulations 9 



12 r- 


P C V*) 


^ /-A 

— PcW + 


4. 8,3 10 r Ia (t) 


14 tvt 

IN 


„£o£ /J.\ 

Pjv w 


= Piv(i) + 


1.1b 1U rj a (t) 


16 n 


Po w 


= PO (*) + 


0.14d ri a (t) 


IVLg 


PMgW 


— PUgW 


+ 8.5 1U J r Ia (t) 


28 Si 


P§?(*) 


= PSi{t) J t 


0.154 r Io (t) 


32g 


Ps'W 


= Ps(t) + 


8.4610- 2 r Io (t) 


40 Ca 


PCa(*) 


= PCa(t) ~ 


f 1.1910- 2 r Io (t) 


56p e 




= PFe(t)~ 


-0.626 r Ja (t) 


z 


p' z ot w 


= Pz(t) + 


1.4 r /a (t) 



Table 4. Total rate of release of chemical yields from a SSP, 
including the contribution of SN la. r i a (t) is given in Table ^ 



Z t (At) = 



Am z 



f: +At e(f)df 



(9) 



4-3.2 Metal release from the star particle 

A particle which has experienced SF is considered to host a 
SSP, which releases gas and metals in time as calculated 
above. From the chemical point of view, that is, within 
such a particle there are "hidden" stellar and gaseous sub- 
components evolving in time. From the point of view of hy- 
drodynamics, however, in our simulations we do not resolve 
these sub-components: a baryonic particle is labelled either 
as gas or as star, and behaves accordingly either as colli- 
sional or collisionless matter. We must hence translate the 
chemical enrichment of a SSP calculated above into chemical 
enrichment of the gas particles. 

Just as in the case of gas restitution and of the SN rates 
(§|4.l| and §4.2), we resort to a statistical approach. A particle 
which has experienced SF may either remain a star forever, 
or become at some point a gas-again particle. As long as it 
remains a star particle, it is assigned the metallicity Zo of 
its SSP; when at some timestep At it turns into a gas-again 
particle, it is assigned the composition Zt(At) of the gas 
released by the SSP within At, given by 

Once more, the main limit of this approach may reside 
in poor statistics: the metal production of a whole stellar 
population is sampled by the, possibly small, fraction of it 
which dies within the timestep At. A large enough number of 
particles is necessary to obtain a meaningful statistical sam- 
pling, and the "chemical timestep" must be chosen appro- 
priat ely, ge nerally longer that the hydro-dynamical timestep 
(see §4.2.3). Our test for the single burst case and our simu- 
lations below indicate that our algorithm works already with 
a few thousand baryonic particles. 



4.4 Metal diffusion in the gas 

Finally, we must describe the chemical evolution of the over- 
all population of gas particles, both gas-again particles and 
those which have never been stars, which proceeds by metal 
diffusion. Particles which were never stars can acquire met- 
als from neighbouring gas-again particles through diffusion. 



Conversely, a gas-again particle is initially assigned the 
metallicity Zt(At) from Eq. ^ as discussed above; from 
then on, its composition will evolve only by metal exchange 
with other nearby gas particles, with no further memory 
that it had been a star in the past. For the sake of metal 
diffusion, that is, once the initial metallicity of a gas-again 
particle has been assigned there is no need to distinguish 
between gas and gas-again particles. 

At each time step, we diffuse metallicity among gas 
particles according to the scheme originally developed by 
Groom (1997), which consists of an SPH translation of the 
usual diffusion equation: 

§— do) 

where k is the diffusion coefficient. In principle, this scheme 
is more physically grounded than the widely used SPH- 
smoothing, in which metals are SPH-spread among the 
neighbouring particles. It relies on the idea that the diffusion 
of metals is driven by SN explosion and energy injection in 
the interstellar medium during the SN remnant phase. 

We derive the diffusion coefficient from the models by 
Thornton et al. (1998). The typical size of a SN remnant 
after 10 6 yr from the explosion is 50-100 pc, whereas the 
typical velocity of the gas accelerated by the SN remnant at 
that time is 40-60 km/sec. In this way we obtain: 

ft = (50 km sec~ 1 )(1.85 x 10 15 km) km 2 sec -1 

The SPH translation, which requires second derivatives of 
the kernel, is 



dZi \ -> 



3=1 



Pi + pi 



(Zi -Zj)x 



(11) 



This SPH algorithm was successfully tested by Groom 
(1997) versus suitable analytical exact counterparts. 

A possible drawback of the adopted diffusion scheme 
could be the use of the 2 n derivative of the smoothing 
function, which may introduce spurious fluctuations in the 
resulting metallicity. The use of a spline kernel, which is con- 
tinuous with its 1 st and 2 n order derivatives, ensures that 
unphysical fluctuations in the diffused quantities (metallic- 
ity in our case) do not appear. 



5 A SINGLE BURST 

A first useful test for our algorithm is to model the evolution 
of the various chemical quantities (gas fraction, metallicity 
etc.) in the simple case of a single burst of SF. This test 
focuses entirely on the implemented chemical network, un- 
coupled in this case from hydrodynamical effects. It allows 
us on one hand to visualize the expected chemical produc- 
tion of a SSP, and on the other hand to check the validity 
of the statistical approach. 

The single burst model is realized with 5000 particles 
of gas, of total mass of 10 11 Mq, turned into stars instanta- 
neously at the beginning of the simulation (To = 0, Zo = 0). 
These star particles are let evolve afterwards according to 
the statistical prescriptions for gas and metal restitution 
given above, without any further SF episodes. 
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Figure 5. Chemical evolution for a single burst of star formation involving 10 11 Mq. Upper-left panel: cumulative rates of SN II 
and la (T < and > 7.5e7 yr, respectively) in logarithm of the number of events per century. Upper-right panel: evolution of the total 
metallicity. Lower panels: evolution of iron (on the left) and oxygen (on the right) abundances. Solid line: numerical results for the 
Salpeter IMF; short-dashedline: Kroupa IMF; long-dashed line: Arimoto-Yoshii IMF. The dotted lines are the corresponding exact 
analytical predictions. 
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Figure 4. Returned gas fraction in time for a single burst of star 
formation involving 5000 particles. Thin dotted lines represent 
the exact analytical expectations. 



It is sufficient to model the evolution of a couple of char- 
acteristic elements with different production timescales, to 
assess the capability of the statistical approach in describing 
detailed chemical evolution in time. We choose to follow the 
evolution of oxygen (representative of a-elements produced 



by massive stars over short time scales), of iron (represen- 
tative of elements with a delayed contribution), and of the 
global metallicity. 

As everything in our algorithm (SN rates, metal pro- 
duction etc.) is entirely governed by the statistical gas resti- 
tution, it is fundamental to check this quantity first of all. 
Fig. ^ shows the evolution of the returned gas fraction in 
time for our single burst test, for the three IMF cases, as 
compared to the exact expectations computed directly from 
the analytical fitting functions. The response of the algo- 
rithm is very good already with ~5000 particles. 

In Fig. ^ we display the results for SN rates and for the 
evolution of metallicity, iron and oxygen abundances in the 
returned gas, again compared to the analytical predictions. 

The results reproduce the expected trends for a SSP, 
and although statistical fluctuations are apparent, for in- 
stance, in the predicted SN rates, the overall behaviour is 
traced quite well. Metallicity is very high at the beginning, 
when the highly metal enriched gas from massive stars is 
expelled. Later, this extremely metal rich gas is diluted by 
more metal poor gas ejected by long-lived stars, so that 
metallicity decreases. The oxygen abundance follows the 
same trend, and represents roughly a half of the global 
metallicity, being oxygen in fact the most abundant of all 
metals. The iron abundance also decreases after the initial 
peak due to massive stars, but contrary to oxygen it later 
stabilizes and even increases again thanks to the contribu- 
tion of SN la. The oxygen-to-iron ratio correspondingly 
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changes from super-solar values at early times, when SN 
II dominate the chemical enrichment, to solar and sub-solar 
values at later times. 

The absolute number of SN is largest for the Arimoto- 
Yoshii IMF, which is the most top-heavy, and decreases 
moving to the Salpeter and to the Kroupa case; the same 
trend is seen in the respective metallicities and oxygen abun- 
dances, dominated by the production of massive stars. The 
relative number of SN la to SN II, however, decreases when 
going from Arimoto- Yoshii to Salpeter to Kroupa because 
these latter IMFs are more skewed toward smaller stellar 
masses. Hence, the final [O/Fe] ratio typical of the Arimoto- 
Yoshii SSP is larger than for the other IMF cases. 

Within the range of massive stars (T < 10 s yr) , one can 
notice that the typical iron abundance in the ejected gas in- 
creases going from the Arimoto- Yoshii to the Salpeter to the 
Kroupa IMF. This is due to the fact that SN from progeni- 
tors in the low mass end of massive stars (say, 8-15 Mq) are 
expected to produce a higher amount of iron, relative to the 
global mass ejection or to the oxygen production, than stars 
in the higher mass range. Hence, the Arimoto- Yoshii IMF 
favours the highest masses with very oxygen-rich, and less 
iron-rich, ejecta; on the contrary, the Kroupa IMF is more 
skewed toward SN with a higher relative iron production. As 
a consequence, the typical [O/Fe] ratio of SN II ejecta from 
an Arimoto- Yoshii SSP is higher than from a Salpeter SSP 
than for a Kroupa SSP. Typical values are [O/Fe]=+0.56 
for Arimoto- Yoshii, +0.50 for Salpeter, +0.45 for Kroupa. 
The theoretical estimates of iron production as a function 
of mass, on the other hand, are very uncertain and this de- 
tailed mass effect might be somewhat spurious. However, all 
of these values are compatible with empirical evidence for 
the a-enhanced stellar population of the Galactic halo (see 
e.g. the data by Carretta et al. 2000 in Fig. p"l| ). 

The Kroupa IMF has been adopted for the simulations 
presented in the following. 



6 A DISC-LIKE GALAXY 

As a first astrophysical test-application of our code we run 
a simulation for an individual galaxy. 

The initial configuration for this object is a a spherical 
DM halo with a density profile oc 1/r (see Lia et al. 2000 for 
a justification of this choice). The total mass of the system 
is 1.3 10 12 Mq with a baryonic fraction equal to 0.1. The 
galaxy is modeled using 30,000 baryonic particles and 15,000 
DM particles. Accordingly, the mass of a DM particle is 
7.8 xlO 7 Af , while the mass of a gas particle is 4.3 x lO 6 Af0 
We assign to the halo a rigid rotation with an angular speed 
A = 0.08. 

Due to both cooling and angular momentum, gas cools 
down in a thin rotating disk forming stars. The distribution 
of star and gas particles are shown in Fig. ^. The snapshots 
refer to 7 Gyr from the beginning of the simulations; while 
star particles are settled into a stellar disc-like object, most 
of the gas is in the form of a hot spheroidal halo. 

In Fig. | we plot the SF history (SFH) of the overall 
object (solid line). The SFR exhibits a strong peak in the 
initial 0.5 Gyrs, then sharply declines and by the age of 
4 Gyrs settles to a rather low constant value of a few Mq /yr. 
Correspondingly, in the final Gyrs of the simulation most of 
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Figure 6. Distribution in the x-y and x-z planes of star particles 
(top panels) and of gas particles (bottom panels) in the innermost 
40 kpc, after 7 Gyrs. 
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Figure 7. Global star formation history of our disc— like galaxy. 
The thin dotted and dashed lines correspond to the low resolution 
and very high resolution simulations, respectively. 



the gas is distributed in a rather hot halo (Fig. |6|) around 
the stellar disc, while little cold gas is left to fuel more active 
SF within the disc. 

To test the behaviour of the "chemical algorithm" in the 
simulation, we compare it with the results of a detailed one- 
zone chemical evolution model with the same SFH: global 
properties like the overall gas consumption, SN rates etc. 
should correspond. Usually, models for chemical evolution 
calculate their own SFH internally after some prescribed 
analytical Schmidt-like law (e.g. Portinari & Chiosi 1999). 
For the present purpose, we developed instead a version of 
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Figure 8. Evolution of the global gas fraction in our disc-like 
galaxy (solid line) compared to the predictions of a detailed chem- 
ical model (dotted line). 
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Figure 9. Evolution of the SN rates. Solid line: type II SN; dashed 
line: type la SN (rate amplified by a factor of 10 for the sake of 
clarity); dotted lines: corresponding predictions from the chemical 
model. 



the chemical model by PCB98 suitable to be force-fed the 
SFH as an input information. The chemical model has then 
been run for a closed system, just as our hydro-dynamical 
galaxy is when considered globally. 

Fig. ^| shows the evolution of the global gas fraction in 
our object, as SF proceeds. The dotted line is the corre- 
sponding prediction from the chemical model, once the SFH 
in Fig. [?] is imposed. The agreement is quite good. 

Fig. ^ shows the time evolution of the SN rates for 
type II and type la SN, respectively. Notice how the rate 
of SN II closely traces the SFR, as expected, while the rate 
of SN la is much more diluted in time. In this object with 
a low final SFR, the two types of SN end up with compa- 
rable rates. The trend for the SN rates, just as for the gas 
consumption in the previous figure, is in full agreement with 
the predictions of the chemical model; this is a confirmation 
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Figure 10. Global (gas+stars) average metal enrichment history 
(dashed lines) and for gas and stars separately (solid and dot- 
ted lines, respectively). Thin lines: SPH simulation; thick lines: 
one-zone chemical model. The global average metal enrichment 
compares very well between the two models (see text). 



of the good behaviour of our algorithm for gas restitution 
and chemical evolution in more complex simulations than a 
single burst. 

With respect to the detailed metallicity evolution of the 
gaseous and stellar components of the system, instead, it is 
little meaningful to compare the hydrodynamical simulation 
and the one-zone chemical model. In fact, a basic assump- 
tion of one-zone models is that the system is always homoge- 
neous in space, namely the metals produced and ejected by 
stars are instantaneously mixed and spread over the whole 
gas mass in the system. This is fundamentally different from 
the behaviour of dynamical models of galaxies, where star 
formation and metal pollution are localized. In a realistic 
disc-like galaxy stars pollute, and form from, the gas in the 
disc region, while there will be a fraction of gas away from 
the disc which may take little or no part to SF and metal 
pollution. Qualitatively, in a dynamical simulation of a disc 
galaxy we expect stars to enrich a more limited amount of 
gas (that located in their surroundings) , which gets enriched 
more effectively than in the one-zone model; this highly en- 
riched gas is somehow "compensated" by the the pristine, or 
almost pristine, gas residing far away from the disc. Conse- 
quently, the stars in the disc form in an environment which 
is more metal enriched than what a simple one-zone model 
predicts, and their metallicity distribution is expected to be 
different than that derived from the chemical model. How- 
ever, if the distribution of metals between gas and stars is 
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Figure 11. Evolution of the average chemical composition of the 
stellar component: [O/Fe] ratio versus metallicity 
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Figure 12. Chemical composition of stars in three different age 
bins (old, intermediate and young). To is the birth time of stars, 
the overall system is 7 Gyr old. 



expected to be different in the dynamical simulation and 
in the chemical model, the global metallicity of the sum of 
the gaseous+stellar components should be comparable be- 
tween the two models, as this traces t he g lobal metal pro- 
duction for that particular SFH. Fig. |l0| illustrates these 
expected effects: the histories of metal enrichment for gas 
and stars separately are evidently very different between 
the one-zone model and the simulation. In particular, in 
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Figure 13. Metallicity profiles for gas (left panels) and stars 
(right panels) 



the hydro-dynamical simulation more metals result locked 
into the stellar component, since the stars form in local- 
ized, very metal enriched regions; correspondingly, the gas 
metallicity in the overall is lower than that of stars because 
the gas component includes also plenty of un-processed hot 
gas in the outskirts of the galaxy. However, when one con- 
siders the global metal enrichment (dashed lines, average 
of gas+stars), which is the record of the overall metal pro- 
duction, the two models compare very well; minor differ- 
ences can be imputed to the fact that the chemical model 
by PCB98 includes detailed metallicity-dependent yields, 
while our fitting formula in Table ^ are based on average 
values of the yields. The comparison in Fig. ^ implies that 
the chemical algorithm in the SPH simulation gives an ade- 
quate description of the metal production. 

Fig. shows the evolution of the global, average [O /Fc] 
ratio for the stellar component, versus metallicity. At low 
metallicities stars display an a-enhanced composition typ- 
ical of SN II enrichment, while at higher metallicities the 
[O /Fe] ratio decreases, as expected from the additional con- 
tribution of SN la to the iron production. Notice how the 
predicted evolution compares to the observational data for 
stars in the Solar Neighbourhood. We remark, however, 
that this comparison must be regarded only as qualitative, 
since the present simulation is not aimed at reproducing 
the Milky Way or the Solar Neighbourhood, but only to 
test the self-consistency of the chemical algorithm we imple- 
mented. Actually, our object in the end resembles more the 
gas-consumed disc of a young SO than that of an Sb spiral 
like our own Galaxy. 
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Fig. |l2j shows the chemical composition of stars in three 
representative age bins: "old" stars born within the first Gyr 
of the simulation, "intermediate age" stars, born between 1 
and 5 Gyrs of the simulation, and "young" stars, i.e. younger 
than 2 Gyr by the end of the simulation. As expected, the 
bulk of old stars concentrate around high [O/Fe] ratios be- 
cause little iron has yet been contributed by SN la , while 
at decreasing age the bulk of stars moves to higher metallic- 
ities and lower [O/Fe] ratios. The stars of the youngest bin 
end up with super-solar metallicity and slightly under-solar 
[a/Fe] composition. 

Finally, in Fig. [l3| we show the metallicity profiles, for 
gas and stars, in the innermost regions at the final age of 
7 Gyr. In these regions, where the bulk of the galaxy resides, 
the gas is much more metal rich than the stars, as expected 
in general since the gas is an instantaneous picture of the 
chemical enrichment reached at present, while stars are a 
signature of the accumulated chemical history starting from 
the initial, metal-poor early epochs. Besides, in the galaxy 
under examination very little gas is left in the late stages 
within the stellar disc, hence its chemical enrichment pro- 
ceeds very fast. 

Metallicity gradients are visible, and they are more pro- 
nounced for the gas component than for the stars, as ex- 
pected in general from chemical evolution models (e.g. Ed- 
munds & Greenhow 1995). 

6.1 Testing the effects of resolution 

A basic issue of hydrodynamical simulations of galaxy for- 
mation is to what extent the results are affected by resolu- 
tion. In principle, this problem might be even more crucial 
when a statistical algorithm is used for SF and chemical 
evolution. 

The resolution effect on the SFH of individual galaxies 
has been explored in Lia et al. (2000) by means of ad hoc 
simulations of the same individual galaxy at increasing num- 
ber of particles (2,000, 20,000 and 200,000). They showed 
that the SF recipe converged above 20,000 particles; beyond 
that, the SFH did not change significantly by varying the 
particle number further. 

We repeat an analogous test here for our disc-like 
galaxy, since our SF algorithm is different now: although 
we adopt the same formal SF law as in Lia et al. (2000), 
this law is now implemented with a probabilistic approach 
(§ |^). Besides, resolution tests are needed not only for the 
convergence of the SFH, but also because the gas and metal 
restitution are treated statistically. 

In the previous section, we discussed the self- 
consistency of our chemical algorithm for a simulation with 
30,000 baryonic particles, by comparing its results with a 
one-zone chemical model, where possible. We will refer to 
this simulation as the "high-resolution" one. Fig. [?] also 
shows the SFH of the same galaxy modelled with "low res- 
olution" (8,000 particles, thin dotted line) and with "very 
high resolution" (200,000 particles, thin dashed line). 

The low resolution case shows a much more noisy SFH, 
although qualitatively the trend resembles that of the high 
resolution case: an initial high SFR is mantained for ~1 Gyr, 
after that its level declines rapidly. 

Fig. [m] illustrates the performance of the chemical al- 
gorithm in the low resolutions case. The gas restitution (top 
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Figure 14. Comparison between the low resolution galaxy simu- 
lation and the chemical one-zone model with the same SFH (the 
dashed line in Fig. |?|). Top panel: evolution of the gas fraction; 
solid and dotted lines for simulation and chemical model respec- 
tively, as in Fig. |i[ Mid panel: metallicity evolution for gas, stars 
and global average; thin and thick lines for simulation and chem- 
ical model respectively, as in Fig. Bottom panel: evolution of 
the SN rates; the dotted lines are the predictions of the chemical 
model. 



panel), and consequently the metal release (mid panel), are 
a bit underestimated when compared to the predictions of 
the one-zone model. This can be imputed to the large oscil- 
lations of the SFR; these oscillations produce even stronger 
noise in the gas (and metal and SN) release, as these imply a 
further "probabilistic event" on top of the one that induces 
SF. This noise is particulary evident in the SN rates (bottom 
panel). However, the effect is not large: gas and metal resti- 
tution seems to be underestimated just by 10% or so. Also 
the SN rates in the simulation, in spite of heavy oscillations, 
on average follow the exact one-zone counterpart. 

The very high resolution simulation was followed up to 
1.75 Gyr, with the sole purpose of demonstrating that the 
SFH converges when the number of particles is increased 
beyond 30,000. As the performances of the chemical algo- 
rithm are very good in the high resolution simulation (see 
previous section), once the convergency of the SFH is also 
demonstrated one can be confident about the overall consis- 
tency of the results. 
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Figure 15. Evolution of the global gas fraction in the cluster 

In summary, concerning resolution effects we can con- 
clude the following. 

• With 8,000 particles the SFH and chemical evolution 
of the object are quite noisy, and the statistical algorithm 
somewhat underestimates the gas and metal release. How- 
ever, for that SFH the mismatch with respect to the detailed 
chemical model is by no means dramatic (within 10%) so the 
gross features of chemical evolution are already rendered. 

• With 30,000 particles and beyond, the SFH converges 
and at that level of resolution the statistical chemical algo- 
rithm also responds pretty well. 

Concerning this latter point, a similar result was found by 
Lia et al. (2000): the SFH in their simulations converged 
beyond 20,000 particles. Since in their case the SF algorithm 
was not statistical, we conclude that resolution effects on the 
SF law itself are dominant, while the probabilistic approach 
has a minor impact on the resolution limit. Most important, 
we stress once more that beyond the resolution limit for 
the SFH, also the statistical chemical algorithm for gas and 
metal restitution performs very well. 



7 A CLUSTER OF GALAXIES 

Other interesting astrophysical applications of chemo- 
hydrodynamical codes lie in cosmological simulations of the 
formation and evolution of clusters, to address the problem 
of the chemical enrichment of the ICM self-consistently. For 
the sake of example, we present here one such simulation. 

The initial conditions were realized by perturbing a 
cubic grid of particles with the displacements field made 



available by the Cluster Com parison Project ( uttp://star- 
www.dur.ac.uk/csf/clusdata/). The initial fluctuation spec- 
trum was taken to have an asymptotic spectral index, n = 1, 
and shape parameter, F = 0.25; the cosmological parame- 
ters assumed were: mean density, £1=1; Hubble constant, 
Ho = 50 kins' 1 Mpc -1 ; present-day linear rms mass fluc- 
tuations in spherical top hat spheres of radius 16 Mpc, 
as = 0.9; and baryon density (in unit of the critical den- 
sity), ttb ~ 0.1. The perturbation was centered on a cu- 
bic region of size L = 64 Mpc. See Frenk et al. (1999) for 
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Figure 16. Chemical enrichment in time for the ICM 



more details. The system was divided in two zones, an inner 
sphere of radius 22 Mpc which was filled with 36 3 gas parti- 
cles and 36 3 dark particles, surrounded by a sphere of radius 
32Afpc in which only dark matter was present. Initially gas 
and dark matter were placed on top of each other and were 
given the same velocities, computed using the Zel'dovich ap- 
proximation and adding the Hubble flow. In the inner region 
the masses of a DM and gas particle are 2.1 x 1O 1O M0 and 
2.4 x 10 9 M Q , respectively. 

Fig. ^ shows the evolution of the gas fraction over the 
total baryonic matter. At the end of the simulation, ~75% 
of the baryons are in gaseous form, while 25% is locked into 
stars, compatible with observed estimates (the mass in gas 
in clusters is 2-5 times that in galaxies, Arnaud et al. 1992). 

Fig. [l(j shows the evolution of chemical abundances in 
the gaseous component. The final metallicity in the ICM is 
~0.25 solar, in broad agreement with observational values. 
Notice how the metalicity evolution in the gas is negligible 
at low redshifts (z < 1), consistent with observational data 
(Mushotzky & Loewenstein 1997, Matsumoto et al. 2000). 
Also the average [O/Fe] ratio is compatible with observa- 
tional data, although these are highly uncertain (Ishimaru 
& Arimoto 1997). 

Fig. [l^ shows the final metallicity distribution in the 
ICM. A hint of the existence of a metallicity gradient (a re- 
cent "hot" issue, De Grandi & Molendi 2001, Finoguenov 
et al. 2000, White 2000) can be seen in Z and in [Fe/H]. 
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Figure 17. Mctallicity gradients for gas in the ICM 
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Figure 19. Radial temperature profile for the , 
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Figure 18. Iron distribution in the ICM: the darkest regions 
indicate zones with higher iron content. Sec fig.jn] for comparison 



In particular, a metallicity peak can be distinguished in the 
very central regions (200-300 kpc) , while further out metal- 
licity seems to be more uniformly distributed, though with a 
large scatter. Concerning the radial behaviour of the [a/Fe] 
ratio, in our simulation there seems to be no differential gra- 
dients for a-elements and iron, as the [O /Fe] ratio seems to 
be roughly uniform at all radii, with slightly subsolar values, 



Figure 20. Global star formation history per unit volume from 
the whole simulation 



though with a large scatter. The iron distribution is shown 
in Fig. while the temperature profile is plotted in Fig. |l9[ 

Finally, in Fig. ^ we show the estimated average evolu- 
tion of the SFR per Mpc 3 for the global simulation. Notice 
how the peak on the SFR is more intense and located at 
earlier phases (z ~ 5) with respect to what observed in the 
field (e.g. Madau et al. 1996, Steidel et al. 1999). This is in 
line with expectations when looking at areas of larger den- 
sity than average; in fact, in our simulation we are dealing 
with a slice of Universe where the formation of a cluster of 
galaxies takes place. 

Higher resolution simulations are needed to assess this 
problem in the necessary detail, resolving also e.g. the gas 
still bound to individual objects from that actually expelled 
into the ICM, the role of ram pressure stripping over indi- 
vidual galaxies in the central parts of the cluster, and so 
forth. Such higher resolution simulations will be addressed 
in future work. In the context of the present paper, we pre- 
sented this simulation here mostly for ilustration purposes, 
with no claims of conclusiveness. 
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8 SUMMARY AND CONCLUSIONS 

In this paper we developed a new algorithm for detailed cal- 
culation of chemical evolution in SPH codes, conceived so 
as to be implemented with minor computational costs into 
simulations with a very large number of particles, and in 
parallel codes. As computer power and capabilities continue 
to improve, allowing for heavier and heavier simulations, the 
opportunity arises to develop algorithms for the "astrophys- 
ical calculations" (star formation, chemical evolution etc.) 
which provide increasing accuracy the larger the number of 
particles involved. At the same time, as the heaviest com- 
putational effort is typically devoted to the calculation of 
gravity forces, these "astrophysical" algorithms should bear 
as least as possible on the performances, mantaining a low 
relative computational load when the number of particles 
increases. Our statistical algorithm for star formation and 
chemical evolution, presented in this paper, is specifically 
addressed to this purpose. 

After testing the self-consistency of the algorithm, with 
two illustrative applications we showed how our statistical 
algorithm, implemented in a SPH code, is suitable to model 
the evolution and distribution of different chemical elements 
in various contexts. It provides a tool to follow both the 
chemo-hydrodynamical evolution of individual galaxies and, 
in the framework of cosmological simulations, the chemi- 
cal enrichment of the ICM/IGM as well as the global SFR 
and cosmic chemical evolution. Besides, the description of 
star particles as independent SSPs of assigned IMF, age and 
metallicity should make it straightforward to follow spectro- 
photometric evolution as well, in the simulations. 

Our algorithm is meant to be easy to implement into 
any SPH code. 
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